% save('data/experiment-60min.mat', 'inst_count_vec', 'rt_vec');
% load('data/experiment-60min.mat');
% load('data/experiment-40min-1servtype.mat');
classdef test_mpc < handle
    
    properties
        % model=[];
        inst_count_vec=[];
        rt_vec=[];
        m=[];
        m_inst_count_vec=[];
        m_rt_vec=[];
        samplingTime = 1;
        lqr=[];
        x0=[];
        system=[]; 
    end
    
    methods
        
        %  clear all
        %  clear
        %  load('data/experiment-60min.mat');
        %  open policy/test_fake_lqr.m
        %  test_fake_lqr().build_model( inst_count_vec, rt_vec)
        function obj=test_mpc(varargin)
            if nargin==1 
                obj.system=varargin{1};   
            end
        end
        
        function [m_x,d_x]=linearize(obj,x)
            m_x = mean(x);
            % d_x = x-m_x;
            d_x = x-ones(size(x,1),1)*m_x;
        end
        
        function d_x=linearize_around(obj,x,m_x)
            d_x = x-ones( size(x,1),1) * m_x ;
        end
        
        function x=delinearize_around(obj,d_x,m_x)
            x = ones( size(d_x,1),1) * m_x  + d_x;
        end
        
        function [opt_instances,opt_rt]=solve_through_lp(obj,sys,x0,t_,m_inst_count_vec,m_rt_vec, cost_weight_r,f_inst)
            parameters=optimset('linprog');
            % OPTIONS=optimset(parameters,'LargeScale','off','Display','off','simplex','on'); %
            OPTIONS=optimset(parameters,'Display','off'); %
            % t_=10;
            % [r1 ... r10, n1,1 ... n1,10,  n2,1,....,n2,10,....]
            % min f'x
            %     f = [ones(t_,1) * cost_weight_r ; ones(t_,1)];
            prices = prices_amz().get_hourly_price();
            for i=1:t_
                f(i,1)=cost_weight_r;
                for j=1:size(m_inst_count_vec,2)
                    f(i+j*t_,1) = prices(j);
                end
            end
            
            % Ax=<b
            A=[];
            b=[];
            
            % lb=<x=<ub
            % 0<r_i<inf
            for i=1:t_
                lb(i,1)=0-m_rt_vec;
                ub(i,1)=Inf-m_rt_vec;
                %  0<n_i<20
                for j=1:size(m_inst_count_vec,2)
                    % set the lower bound according to the reserved amount f_inst(j,i)
                    lb(i+j*t_,1) =  0-m_inst_count_vec(j); %f_inst(j,i)-m_inst_count_vec(j);  
                    ub(i+j*t_,1)=20-m_inst_count_vec(j);
                end
            end
            
            % Aeq.x=beq
            Aeq=zeros(t_-1 , (size(m_inst_count_vec,2)+1)*t_ );
            beq=zeros(t_-1,1);
            for i=1:t_-1
                % r_i = a.r_(i-1)+b.(n_i);
                % -a.r_(i-1)+1*r_i -b.(n_i)=0;
                Aeq(i,i) = -1 * sys.a;
                Aeq(i,i+1) = 1;
                b_=sys.b;
                for j=1:size(m_inst_count_vec,2)
                    Aeq(i,i+j*t_) = -1 * b_(j);
                end
            end
            % r1=x0;
            Aeq=[Aeq ; [1 zeros(1,((size(m_inst_count_vec,2)+1)*t_-1))]];
            beq=[beq ; x0];
            
            % lb=[]; ub=[];
            
            x0=[];
            % [x,fval,exitflag]      = linprog(f,         A ,b,  Aeq ,beq                     ,lb               ,ub  ,x0  ,options)
            [x,f2,flag2]              = linprog(f, A ,b, Aeq  ,beq , lb , ub    , x0    ,OPTIONS);
            
            
            if flag2 == -3
                disp('unbounded');
                disp(' ');
            else
               % disp([f2,flag2]);
                opt_instances = [reshape(x(t_+1:end,1),t_,size(m_inst_count_vec,2)) + ones(t_,1)*m_inst_count_vec]';
                opt_rt = [x(1:t_,1) + m_rt_vec]';
            end
            
        end
        
        function build_model(obj, inst_count_vec, rt_vec) 
            load('data/experiment-60min.mat'); 
            obj.inst_count_vec=inst_count_vec;
            obj.rt_vec=rt_vec;
            %----------
            n_st=1;% 1 state variable
            n_in= size(inst_count_vec,2);  %11 inputs
            n_out= 1; %size(rt_vec,2); % 1 output
            n_env=1; % environmental disturbances
            A=ones(n_st,n_st);
            B=ones(n_st,n_in);
            K=ones(n_st,n_env);
            C=eye(n_out,n_st);
            D=zeros(n_out,n_in);
            
            
            [m_rt_vec,              d_rt_vec]= obj.linearize(rt_vec);
            [m_inst_count_vec,  d_inst_count_vec] = obj.linearize(inst_count_vec);
            
            
            %             m_rt_vec=0; %2
            %             [  d_rt_vec]= obj.linearize_around(rt_vec , m_rt_vec);  %  2  seconds desired
            %              m_inst_count_vec= ones(1, size(inst_count_vec,2))* 0; %2
            %           [ d_inst_count_vec] = obj.linearize_around(inst_count_vec,...
            %                         m_inst_count_vec);
            
            u =  d_inst_count_vec;
            y = d_rt_vec;
            
            data=iddata(y, u , obj.samplingTime);
            % m=idss(A,B,C,D,K,'Ts', obj.samplingTime,'DisturbanceModel','Estimate','InitialState','Estimate');
            m=idss(A,B,C,D,K,'Ts', obj.samplingTime,'DisturbanceModel','None','InitialState','Estimate');
            %             m0.As = [NaN,0;0,NaN];
            %             m0.Bs = [NaN;NaN];
            
            m.Cs = C;
            m1=init(m);
            mm=pem(data,m1);
            
            obj.m=mm;
            obj.m_inst_count_vec=m_inst_count_vec;
            obj.m_rt_vec = m_rt_vec;
            
            %-------------------------------------------------
        end
           
        
        function [instance_quantity_to_fire,instance_quantity_to_destroy]=get_choice_at_time(obj,t)
            t_=20;
            control_affect_window=5; 
            cost_weight_r=1000000;
            sys = ss(obj.m,'m');
            f_inst = obj.system.future_inst()';
            % grab info on already allocated instances 
            if size(f_inst,2)<t_
                f_inst=[f_inst zeros(11,t_-size(f_inst,2))];           
            end
            % update the model state based on system response 
            past_rt=obj.system.past_util();
            if size(past_rt,2)>0
                obj.m.x0=past_rt(end);    
            end                        
            [opt_instances,opt_rt] = obj.solve_through_lp(sys,obj.m.x0,t_, obj.m_inst_count_vec, obj.m_rt_vec, cost_weight_r,f_inst);       
            instance_quantity_to_destroy=zeros(1,11); 
            % see how much instances we need more than allocated in the
            % control_affect_window. then take the mean over the window 
            % and allocate the extra         
            instance_quantity_to_fire = opt_instances(:,1:control_affect_window)...
                                            - f_inst(:,1:control_affect_window);
   %         instance_quantity_to_fire=abs(instance_quantity_to_fire);
              instance_quantity_to_fire = mean(instance_quantity_to_fire,2);  
             instance_quantity_to_fire(instance_quantity_to_fire<0)=0;
     
            instance_quantity_to_fire = round(instance_quantity_to_fire)'; 
            % see system current number of inst.....
            kk=0;
            % compare 
            
            ......
            
%             default=[0 0 0 0 0 0 0 0 0 0 0];
%             % min, type, num
%             instance_quantity_to_fire=default;
%             instance_quantity_to_destroy=default;
%             t=[0,1];
            %             instance_quantity_to_fire(type)=num;
            %             instance_quantity_to_destroy(type)=num;
            
        end
        
         % ------------------------------------------ test -----------------------------------
            %             load('data/experiment-60min.mat');
            %            load('data/experiment-40min-1servtype.mat');
            % %            noise=randn(60,2)*0.001;
            %             x_sys = mm.X0; % system state  for the first time
            %             rt_acc=[];
            %             u_opt_acc=[];
            %             for i=1:20
            %                 %------ system update ------
            %                 x_sys =sys.a*x_sys + sys.b* (opt_instances(:,i)-m_inst_count_vec') ;
            %                 rt_acc=[rt_acc x_sys+m_rt_vec];
            %                 % u_opt_acc=[u_opt_acc ; u_opt'];
            %             end
            %             rt_acc
            
            %             reference = [ones(1,20)*8 ones(1,20)*6
            %             ones(1,20)*7];
            %             d_reference = obj.linearize_around(reference , m_rt_vec);  %  2  seconds desired
            %
            %
            %             t= 0:1:60-1;
            %
            %          [y,t,x] =    lsim(sys_all, [d_reference ; input] ,t);
            %          figure
            %          u=UtilityLib();
            %          [0.4,0.4,0.4]
            %          % plot(reference, strcat(u.lnstyle{2}),...
            %          plot(d_reference, strcat(u.lnstyle{2}),...
            %                 'LineWidth',2,...
            %                 'MarkerSize',10);
            %             hold on;
            %         % plot(obj.delinearize_around(y(:,1) , m_rt_vec), strcat(u.lnstyle{3},u.color{7}),...
            %          plot(y(:,1) , strcat(u.lnstyle{3},u.color{7}),...
            %                 'LineWidth',2,...
            %                 'MarkerSize',10);
            %          legend('Reference','System')
            %          title('System/Reference response time');
            %
            %          figure;
            %          alloc_instances =  obj.delinearize_around(y(:,2:12) , m_inst_count_vec);
            %          % plot_instances(alloc_instances); legend(prices_amz().get_api_names');
            %          plot_instances(y(:,2:12)); legend(prices_amz().get_api_names');
            %          legend(prices_amz().get_api_names')
            %          title('Number of Instances');
            %
            
    end
    
end

